/*******************************************************************************
*
*   Copyright ©2009 Alexandru Marin
*
*   This file is part of nwave/pMRI C library.
*
*   This project started as an MSc degree project during a
*   stage at the Paris-Est university, under joint supervision of prof.
*   Bogdan Ionescu (Image Analysis and Processing Laboratory, "Politehica"
*   University of Bucharest), Caroline Chaux and Jean-Christophe Pesquet
*   (Laboratoire d'Informatique, Université Paris-Est, Institut Gaspard Monge).
*
*   nwave/pMRI C libray is free software: you can redistribute it and/or modify
*   it under the terms of the GNU General Public License as published by
*   the Free Software Foundation, either version 2 of the License, or
*   (at your option) any later version.
*
*   nwave/pMRI C libray is distributed in the hope that it will be useful,
*   but WITHOUT ANY WARRANTY; without even the implied warranty of
*   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
*   GNU General Public License for more details.
*
*   You should have received a copy of the GNU General Public License
*   along with nwave/pMRI C library.  If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/


#include "mri_utils.h"
#include <math.h>
#include <stdarg.h>
#include <string.h>

static clock_t time_mark = 0;

void mri_tic() {
	time_mark = clock();
	printf("* * * MRI timer start... * * *\n");
	fflush(stdout);
}

void mri_toc() {
	clock_t diff = clock() - time_mark;
	printf("* * * MRI Timer Stop. %f seconds elapsed. * * *\n", (diff / (float) CLOCKS_PER_SEC));
	fflush(stdout);
}

void mri_write_complex(FILE* f, complex float x) {
	assert(f!=NULL);
	fprintf(f, "%f%c%fi", crealf(x), cimagf(x) > 0 ? '+' : '-', fabs(cimagf(x)));
}

void mri_write_complex_matrix(FILE* f, complex float* mat, mri_matrix_descriptor_t* desc) {

	assert(f!=NULL);
	assert(mat!=NULL);

//	printf("* * * mri_write_complex_matrix \n");

	int totalsize = 1, i;

	// Compute total size
	for (i = 0; i < desc->ndims; i++)
		totalsize *= desc->dims[i];

	// Write header. First number of dimensions
	fprintf(f, "%d", desc->ndims);
	fputc(',', f);
	// Then all dimensions
	for (i = 0; i < desc->ndims; i++) {
		fprintf(f, "%d", desc->dims[i]);
		fputc(',', f);
	}

	// Write actual data
	for (i = 0; i < totalsize; i++) {
		mri_write_complex(f, mat[i]);
		if (i != totalsize - 1)
			fputc(',', f);
	}
}

void mri_write_real_matrix(FILE* f, float* mat, mri_matrix_descriptor_t* desc) {

	assert(f!=NULL);
	assert(mat!=NULL);

	int totalsize = 1, i;

	// Compute total size
	for (i = 0; i < desc->ndims; i++)
		totalsize *= desc->dims[i];

	// Write header. First number of dimensions
	fprintf(f, "%d", desc->ndims);
	fputc(',', f);
	// Then all dimensions
	for (i = 0; i < desc->ndims; i++) {
		fprintf(f, "%d", desc->dims[i]);
		fputc(',', f);
	}

	// Write actual data
	for (i = 0; i < totalsize; i++) {
		fprintf(f, "%f", mat[i]);
		if (i != totalsize - 1)
			fputc(',', f);
	}
}

int mri_read_complex(FILE* f, complex float * result) {

	assert(f!=NULL);
	assert(result!=NULL);

	float fl;
	char c;

	if (fscanf(f, "%f", &fl) == EOF)
		return EOF;

	c = fgetc(f);

	if (c == '+' || c == '-') {
		// Was real part of a complex number
		*result = fl+0*I;
		// Turn back one position (sign)
		ungetc(c, f);
		// Read imaginary part
		if (fscanf(f, "%f", &fl) == EOF)
			return EOF;
		*result += fl * I;
		// Read the i
		fgetc(f);
	} else if (c == 'i') {
		// Was an imaginary number
		*result = fl * I;
	} else {
		// Was a real number
		*result = fl+0*I;
		// In case we haven't reached the eof, rewind
		if (c != EOF)
			ungetc(c,f);
	}

	// Success exit
	return 1;
}

complex float* mri_read_complex_matrix(FILE* f, mri_matrix_descriptor_t* desc) {

	assert(f!=NULL);
	assert(desc!=NULL);

	int i, totalsize = 1;

	complex float* result;

	// Read number of dimensions
	if (fscanf(f, "%d,", &(desc->ndims)) != 1)
		return 0;

	// Allocate dimensions and gaps vector
	desc->dims = (int*) malloc(desc->ndims * sizeof(int));
	assert(desc->dims != NULL);

	desc->gaps = (int*) malloc(desc->ndims * sizeof(int));
	assert(desc->gaps != NULL);

	// Read matrix dimensions and compute gaps
	for (i = 0; i < desc->ndims; i++) {
		if (fscanf(f, "%d,", desc->dims + i) != 1)
			return 0;

		totalsize *= desc->dims[i];

		if(i==0)
			desc->gaps[i]=1;
		else
			desc->gaps[i]=desc->dims[i-1]*desc->gaps[i-1];
	}

	// Allocate memory
	result = (complex float*) malloc(totalsize * sizeof(complex float));
	assert(result != NULL);

	// Read matrix
	for (i = 0; i < totalsize; i++) {
		if (mri_read_complex(f, result + i) != 1)
			break;
		// Read in the delimiter
		fgetc(f);
	}

	if (i != totalsize) {
		// We haven't read up all the matrix, so we free the mem and return error
		free(result);
		return NULL;
	}

	// All fine, return success
	return result;
}

float* mri_read_real_matrix(FILE* f, mri_matrix_descriptor_t* desc) {

	assert(f!=NULL);
	assert(desc!=NULL);

	int i, totalsize = 1;

	float* result;

	// Read number of dimensions
	if (fscanf(f, "%d,", &(desc->ndims)) != 1)
		return 0;

	// Allocate dimensions and gaps vector
	desc->dims = (int*) malloc(desc->ndims * sizeof(int));
	assert(desc->dims != NULL);

	desc->gaps = (int*) malloc(desc->ndims * sizeof(int));
	assert(desc->gaps != NULL);

	// Read matrix dimensions and compute gaps
	for (i = 0; i < desc->ndims; i++) {
		if (fscanf(f, "%d,", desc->dims + i) != 1)
			return 0;

		totalsize *= desc->dims[i];

		if(i==0)
			desc->gaps[i]=1;
		else
			desc->gaps[i]=desc->dims[i-1]*desc->gaps[i-1];
	}

	// Allocate memory
	result = (float*) malloc(totalsize * sizeof(float));
	assert(result != NULL);

	// Read matrix
	for (i = 0; i < totalsize; i++)
		// NOTE that "%f," *will* match the last element, even if it is not followed by ','
		if (fscanf(f, "%f,", result + i) != 1)
			break;

	if (i != totalsize) {
		// We haven't read up all the matrix, so we free the mem and return error
		free(result);
		return NULL;
	}

	// All fine, return success
	return result;

}

char* MRI_PATH = NULL;

void mri_set_path(char* path) {

	if (MRI_PATH != NULL)
		free(MRI_PATH);

	MRI_PATH = (char*) malloc((strlen(path) + 1) * sizeof(char));
	assert(MRI_PATH!=NULL);

	strcpy(MRI_PATH, path);
	printf("* MRI: MRI path set to: %s\n", MRI_PATH);

}

char* mri_get_path() {

	if (MRI_PATH == NULL)
		mri_set_path("./");

	return MRI_PATH;

}

/**
 *  Returns a new string containing all given strings concatenated.
 */
char* mri_strcat(int nstr, char* str1, ...) {

	// Compute result size
	int totalsize = strlen(str1);
	int i;
	va_list ap;
	va_start(ap, str1);
	for(i=1;i<nstr;i++)
		totalsize += strlen(va_arg(ap, char*));
	va_end(ap);

	// Allocate
	char* result = (char*)malloc((totalsize + 1) * sizeof(char));

	// Copy parts
	strcpy(result, str1);

	va_start(ap, str1);
	for(i=1;i<nstr;i++)
		strcat(result, va_arg(ap, char*));
	va_end(ap);

	return result;
}

float* mri_load_real_matrix(char* matrixname, mri_matrix_descriptor_t* desc) {

	char* fullpath = mri_strcat(3, mri_get_path(), matrixname, ".csv");

	FILE* f = fopen(fullpath, "r");
	assert(f!=NULL);

	float* result = mri_read_real_matrix(f, desc);

	fclose(f);
	free(fullpath);

	return result;
}

complex float* mri_load_complex_matrix(char* matrixname, mri_matrix_descriptor_t* desc) {

	char* fullpath = mri_strcat(3, mri_get_path(), matrixname, ".csv");

	FILE* f = fopen(fullpath, "r");
	assert(f!=NULL);

	complex float* result = mri_read_complex_matrix(f, desc);

	fclose(f);
	free(fullpath);

	return result;
}

void mri_save_real_matrix(char* matrixname, float* mat, mri_matrix_descriptor_t* desc) {

	char* fullpath = mri_strcat(3, mri_get_path(), matrixname, ".csv");

	FILE* f = fopen(fullpath, "w");
	assert(f!=NULL);

	mri_write_real_matrix(f, mat, desc);

	fclose(f);
	free(fullpath);
}

void mri_save_complex_matrix(char* matrixname, complex float* mat, mri_matrix_descriptor_t* desc) {

	char* fullpath = mri_strcat(3, mri_get_path(), matrixname, ".csv");

	FILE* f = fopen(fullpath, "w");
	assert(f!=NULL);

	mri_write_complex_matrix(f, mat, desc);

	fclose(f);
	free(fullpath);
}



